Quantile regression

Theory, implementations, and applications

Anton Antonov

Accendo Data, LLC

Abstract

In this talk we are going to discuss the theory, implementation, and application of quantile regression. We can say that least squares linear regression corresponds to finding the mean of a single distribution. Similarly, quantile regression corresponds to finding quantiles of a single distribution. Quantile regression is more robust than linear regression and using quantile regression curves we can obtain more complete pictures of the distributions we analyze.

In the talk we are going to discuss the theoretical background of quantile regression computation, and several implementations: with user specified functions and with B-splines, through linear programming and global minimization. Several applications will be shown for weather forecast, non-linear data fitting, time series analysis, and time series outlier detection. Extensions and future directions are will be discussed.

Examples

Here is an example of quantile and linear regression through sinusoidal noise over a parabola:

Here is a link to a movie:

Time series conversational engine: http://www.youtube.com/watch?v=wlZ5ANglVI4 ,

that demonstrates a conversational engine for time series analysis that utilizes quantile regression.

References

[1] Anton Antonov, Quantile regression Mathematica package, source code at GitHub, https://github.com/antononcube/MathematicaForPrediction, package QuantileRegression.m, (2013).

[2] Roger Koenker, Gilbert Bassett Jr., “Regression Quantiles”, Econometrica, 46(1), 1978, pp. 33-50. 
JSTOR URL: http://links.jstor.org/sici?sici=0012-9682%28197801 %2946 %3 A1 %3 C33 %3 ARQ %3 E2 .0.CO%3 B2-J .

[3] Wikipedia, Quantile regression, http://en.wikipedia.org/wiki/Quantile_regression .

[4] Brian Cade, Barry  Noon, “A gentle introduction to quantile regression for ecologists”, Front. Ecol. Environ. 1(8), 2003, pp. 412–420.

[5] Roger Koenker, Quantile regression, ‪Cambridge University Press‬ - ‪Business & Economics, ‪May 9, 2005.‬

Some of my work

In the past

Kernel developer at Wolfram Research, Inc. (Numerical algorithms)

Principal engineer at Sezmi, Inc. ( TiVo ⊗ Netflix )

Senior research scientist at 
- Panasonic Automotive Systems America (Route prediction, recommendations, conversational engines)
- KIT Digital (aka Piksel)
- Elance

Operations research scientist at RaceTrac Petroleum (logistics optimization and consumption prediction in retail)

Now

Senior research scientist at Senzari, Inc. ( music recommendations, song similarities, playlist analysis, mining, and construction )

Introduction

We can say that least squares linear regression corresponds to finding the mean of a single distribution. Similarly, quantile regression corresponds to finding quantiles of a single distribution. With quantile regression we obtain curves — “regression quantiles” — that together with the least squares regression curve would give a more complete picture of the distributions (the Y’s) corresponding to a set of X’s.

The Quantile Regression Problem (QRP) is formulated as a minimization of the sum of absolute differences. Further, QRP is re-formulated as a linear programming problem, which allows for efficient computation. To compute quantiles other than the median the so called “tilted” function is used. (See this Wikipedia entry.)

For a complete, interesting, and colorful introduction and justification to quantile regression I refer to this article:

Roger Koenker, Gilbert Bassett Jr., "Regression Quantiles", Econometrica, Vol. 46, No. 1. (Jan., 1978), pp. 33-50. JSTOR URL: http://links.jstor.org/sici?sici=0012-9682%28197801%2946%3A1%3C33%3ARQ%3E2.0.CO%3B2-J .

I recently implemented a Mathematica package for computing regression quantiles using Mathematica's function LinearProgramming — see the package QuantileRegression.m provided by the MathematicaForPrediction project at GitHub. Also see this guide for using the function QuantileRegressionFit provided by the package. The guide has some theoretical explanations and shows how the quantile regression problem can be formulated as a linear programming problem.

Conversational engine for time series analysis

http://www.youtube.com/watch?v=wlZ5ANglVI4

Theoretical background

We can formulate the quantile regression problem in way analogous to the formulation of least squares (conditional mean) regression.

Consider a random variable Y having some distribution function F and a sample (Subscript[{(Subscript[y, i])}, i=1]^n) of Y. The median of the set of samples (Subscript[{(Subscript[y, i])}, i=1]^n) can be defined as the solution of the minimization problem

(min_β) (∑_(i=1)^n)|(Subscript[y, i])-β|, β∈ℝ .

To see that the β which minimizes (1) is the median of (Subscript[{(Subscript[y, i])}, i=1]^n), consider two points (Subscript[y, 1])<(Subscript[y, 2]). Then |(Subscript[y, 1])-μ|+|μ-(Subscript[y, 2])|=|(Subscript[y, 1])-(Subscript[y, 2])|, ∀μ∈[(Subscript[y, 1]),(Subscript[y, 2])], hence any μ∈[(Subscript[y, 1]),(Subscript[y, 2])] minimizes (1). Using the observation for two points we see that for three points (Subscript[y, 1])<(Subscript[y, 2])<(Subscript[y, 3]), μ=(Subscript[y, 2]) minimizes (1). For four points (Subscript[y, 1])<(Subscript[y, 2])<(Subscript[y, 3])<(Subscript[y, 4]) any μ∈[(Subscript[y, 2]),(Subscript[y, 3])] minimizes (1). We can generalize these observations and show that (1) gives the median for any set of points.

If we want to find θ-th sample quantile of (Subscript[{(Subscript[y, i])}, i=1]^n) then we need to change (1) into

(min_β)( (∑_((Subscript[y, i])≥β))θ|(Subscript[y, i])-β|+(∑_((Subscript[y, i])<β))(1-θ)|(Subscript[y, i])-β| ), β∈ℝ.

Consider a set of random variables (Subscript[Y, i]), i∈[1,n], n∈ℕ that are paired with a set of x-coordinates X=(Subscript[{(Subscript[x, i])}, i=1]^n). We have data of pairs (Subscript[{(Subscript[x, i]),(Subscript[y, i])}, i=1]^n), where (Subscript[y, i]) is a realization of (Subscript[Y, i]).

The linear regression problem can be formulated as

(min_((Subscript[β, 0]),(Subscript[β, 1]))) (∑_(i=1)^n)((((Subscript[y, i])-((Subscript[β, 0])+(Subscript[β, 1]) (Subscript[x, i]))))^2).

Similarly, the median regression problem can be formulated as

(min_((Subscript[β, 0]),(Subscript[β, 1]))) (∑_(i=1)^n)|(Subscript[y, i])-((Subscript[β, 0])+(Subscript[β, 1]) (Subscript[x, i]))| .

and the θ-th quantile regression problem as

(min_((Subscript[β, 0]),(Subscript[β, 1])))( (∑_(i∈{i:(Subscript[y, i])≥(Subscript[β, 0])+(Subscript[β, 1])(Subscript[x, i])}))θ|(Subscript[y, i])-((Subscript[β, 0])+(Subscript[β, 1])(Subscript[x, i]))|+(∑_(i∈{i:(Subscript[y, i])<(Subscript[β, 0])+(Subscript[β, 1])(Subscript[x, i])}))(1-θ)|(Subscript[y, i])-((Subscript[β, 0])+(Subscript[β, 1])(Subscript[x, i]))| ), (Subscript[β, 0]),(Subscript[β, 1])∈ℝ.

In order to convert (5) into a linear programming problem, let us introduce the non-negative variables (Subscript[u, i]) and (Subscript[v, i]) for which the following equations are true:

(Subscript[y, i])-((Subscript[β, 0])+(Subscript[β, 1]) (Subscript[x, i]))+(Subscript[u, i])=0, i∈{i:(Subscript[y, i])≥(Subscript[β, 0])+(Subscript[β, 1])(Subscript[x, i])},
(Subscript[u, i])=0, i∉{i:(Subscript[y, i])≥(Subscript[β, 0])+(Subscript[β, 1])(Subscript[x, i])},

((Subscript[β, 0])+(Subscript[β, 1]) (Subscript[x, i]))-(Subscript[y, i])+(Subscript[v, i])=0, i∈{i:(Subscript[y, i])<(Subscript[β, 0])+(Subscript[β, 1])(Subscript[x, i])},
(Subscript[v, i])=0,i∉{i:(Subscript[y, i])<(Subscript[β, 0])+(Subscript[β, 1])(Subscript[x, i])}.

Since (Subscript[u, i]) and (Subscript[v, i]) are greater than 0 on complementary sets, we can re-write (6) and (7) simply as

(Subscript[y, i])-((Subscript[β, 0])+(Subscript[β, 1]) (Subscript[x, i]))+(Subscript[u, i])-(Subscript[v, i])=0,(Subscript[u, i])≥0,(Subscript[v, i])≥0,i∈[1,n].

Then (5) expressed with (Subscript[u, i]) and (Subscript[v, i]) becomes

(min_((Subscript[u, i]),(Subscript[v, i]),(Subscript[β, 0]),(Subscript[β, 1])))( (∑_(i∈{i:(Subscript[y, i])≥(Subscript[β, 0])+(Subscript[β, 1])(Subscript[x, i])}))θ (Subscript[u, i])+(∑_(i∈{i:(Subscript[y, i])<(Subscript[β, 0])+(Subscript[β, 1])(Subscript[x, i])}))(1-θ) (Subscript[v, i])) .

Since (Subscript[u, i])≥0 and (Subscript[v, i])≥0 the minimization function (9) can simply be written as

(min_((Subscript[u, i]),(Subscript[v, i]),(Subscript[β, 0]),(Subscript[β, 1])))( (∑_(i=1)^n)θ (Subscript[u, i])+(∑_(i=1)^n)(1-θ) (Subscript[v, i])) .

The equations (8) and formula (10) are the linear programming formulation of the quantile regression problem (5).

Note that (Subscript[u, i])(Subscript[v, i])=0, ∀i∈[1,n].

The quantile regression formulations (5), and (8) and (10) can be done for any model of (Subscript[Y, i]) that is a linear combination of functions over X not just for the linear model (Subscript[β, 0])+(Subscript[β, 1])X.

QuantileRegressionFit

Logarithmic curve with noise

Let us generate some data.

data=LogarithmicCurveWithNoise[1200,10.,100.];
ListPlot[data,AxesLabel->{"x","y"},PlotRange->All,ImageSize->400]
AutoCollapse[]

Consider the following quantiles:

qs={0.05,0.25,0.5,0.75,0.95};

We want to find curves that separate the data according the quantiles. Those curves are called “regression quantiles”.

Pretending that we do not know how the data is generated, just by looking at the plot we assume that the model for the data is

y=(Subscript[β, 0])+(Subscript[β, 1]) x+(Subscript[β, 2]) (Sqrt[x])+(Subscript[β, 3]) log(x).

Let us put the model functions for the regression fit in the variable funcs:

funcs={1,x,(Sqrt[x]),Log[x]};

Here we find the regression quantiles:

qrFuncs=QuantileRegressionFit[data,funcs,x,qs];
TableForm[List/@qrFuncs]

We also apply Fit to the data and the model functions in order to compare the regression quantiles with the least-squares regression fit:

fFunc=Fit[data,funcs,x]

Here is a plot that combines the found regression quantiles and least squares fit:

grData=ListPlot[data,PlotRange->All];grQReg=Plot[Evaluate[MapThread[Tooltip[#1,#2]&,{Append[qrFuncs,fFunc],Append[Map["ϱ("<>ToString[#]<>",x)"&,qs],"Least squares"]}]],{x,Min[data[[All,1]]],Max[data[[All,1]]]},PlotStyle->Append[Table[{AbsoluteThickness[2],Blend[{Yellow,Darker[Red]},i/Length[qs]]},{i,Length[qs]}],{AbsoluteThickness[2],Darker[Green]}],PlotLegends->SwatchLegend[Append[Map["ϱ("<>ToString[#]<>",x)"&,qs],"Least squares"]]];
Show[{grData,grQReg},ImageSize->400]
AutoCollapse[]

Let us check how good the regression quantiles are for separating the data according to the quantiles they were computed for:

tbl=Table[{qs[[i]],Length[Select[data,#[[2]]≥(qrFuncs[[i]]/.x->#[[1]])&]]/Length[data]//N},{i,Length[qs]}];
TableForm[tbl,TableHeadings->{None,{"quantile","fraction\nabove"}}]
AutoCollapse[]

Data generation definitions

Clear[LogarithmicCurveWithNoise]
LogarithmicCurveWithNoise[nPoints_Integer,start_?NumberQ,end_?NumberQ]:=
Block[{data},
data=Table[{t,5+Log[t]+RandomReal[SkewNormalDistribution[0,Log[t]/5,12]]},{t,Rescale[Range[1,nPoints],{1,nPoints},{start,end}]}];
data
];

Clear[SinWithUpwardTrend]
SinWithUpwardTrend[nPoints_Integer,start_?NumberQ,end_?NumberQ]:=
Block[{data},
data=Table[{t,Sin[t π/15]+RandomReal[SkewNormalDistribution[0,Log[t]/3,12]]},{t,Rescale[Range[1,nPoints],{1,nPoints},{start,end}]}];
data[[All,2]]=data[[All,2]]+data[[All,1]]^(1/2);
data
];

Clear[SinWithParabolaTrend]
SinWithParabolaTrend[nPoints_Integer,start_?NumberQ,end_?NumberQ]:=
Block[{data},
data=Table[{t,Sin[t π/10]+RandomReal[NormalDistribution[0,(Sqrt[t])/10]]},{t,Rescale[Range[1,nPoints],{1,nPoints},{start,end}]}];
data[[All,2]]=data[[All,2]]+(Mean[{start,end}]-data[[All,1]])^2/300;
data
];

Quantile regression robustness

Outliers above a constant

Let us demonstrate the robustness of the regression quantiles with the data of the previous example. Suppose that for some reason 50% of the data y-values greater than 10.8 are altered by multiplying them with a some greater than 1 factor, say, α=1.8 . Then the altered data looks like this:

α=2.8;
dataAlt=Map[If[RandomInteger[{0,1}]==1&&#[[2]]>10.5,{#[[1]],α #[[2]]},#]&,data];
ListPlot[dataAlt,AxesLabel->{"x","y"},PlotRange->All,ImageSize->400]

Let us compute the regression quantiles for the altered data:

qrFuncsAlt=QuantileRegressionFit[dataAlt,funcs,x,qs];
TableForm[List/@Chop[qrFuncs]]

and let us also compute the least squares fit of the model (1):

fFunc=Fit[data,funcs,x]
fFuncAlt=Fit[dataAlt,funcs,x]

Here is a plot of the altered data and all fitted functions:

grData=ListLogPlot[dataAlt,PlotRange->All];grQReg=LogPlot[Evaluate[MapThread[Tooltip[#1,#2]&,{Append[qrFuncsAlt,fFuncAlt],Append[Map["ϱ("<>ToString[#]<>",x)"&,qs],"Least squares"]}]],{x,Min[data[[All,1]]],Max[data[[All,1]]]},
PlotRange->All,PlotStyle->Append[Table[{AbsoluteThickness[2],Blend[{Yellow,Darker[Red]},i/Length[qs]]},{i,Length[qs]}],{AbsoluteThickness[2],Darker[Green]}],PlotLegends->SwatchLegend[Append[Map["ϱ("<>ToString[#]<>",x)"&,qs],"Least squares"]]];
Show[{grData,grQReg},ImageSize->400]
AutoCollapse[]

We can see that the new regression quantiles computed for 0.05, 0.25, and 0.5 have not changed significantly:

Grid[List/@Flatten[Transpose[{Chop[qrFuncs],Chop[qrFuncsAlt]}]],Alignment->Left,Dividers->{{True,True},Flatten@{Table[{True,False},{Length[qrFuncs]}],True}}]
AutoCollapse[]

ant that they are still good for separating the un-altered data:

tbl=Table[{qs[[i]],Length[Select[data,#[[2]]≥(qrFuncsAlt[[i]]/.x->#[[1]])&]]/Length[data]//N},{i,Length[qs]}];TableForm[tbl,TableHeadings->{None,{"quantile","fraction\nabove"}}]
AutoCollapse[]

Also we can see that the least squares fit of (11) has significantly changed:

fFunc

fFuncAlt

Quantile regression robustness

Outliers above a chosen regression quantile

Let us demonstrate the robustness of the regression quantiles with different type of data alternation. Suppose that for some reason 70% of the data y-values above the last quantile function are altered by multiplying them with a some greater than 1 factor, say, α=10 . Then the altered data looks like this:

α=10.;
dataAlt=Map[If[RandomReal[{0,1}]>0.3&&(#[[2]]>(qrFuncs[[5]]/.x->#[[1]])),{#[[1]],α #[[2]]},#]&,data];
ListLogPlot[dataAlt,AxesLabel->{"x","y"},PlotRange->All,ImageSize->400]

Let us compute the regression quantiles for the altered data:

qrFuncsAlt=QuantileRegressionFit[dataAlt,funcs,x,qs];
TableForm[List/@Chop[qrFuncs]]

and let us also compute the least squares fit of the model (1):

fFunc=Fit[data,funcs,x]
fFuncAlt=Fit[dataAlt,funcs,x]

Here is a plot of the altered data and all fitted functions:

grData=ListLogPlot[dataAlt,PlotRange->All];grQReg=LogPlot[Evaluate[MapThread[Tooltip[#1,#2]&,{Append[qrFuncsAlt,fFuncAlt],Append[Map["ϱ("<>ToString[#]<>",x)"&,qs],"Least squares"]}]],{x,Min[data[[All,1]]],Max[data[[All,1]]]},
PlotRange->All,PlotStyle->Append[Table[{AbsoluteThickness[2],Blend[{Yellow,Darker[Red]},i/Length[qs]]},{i,Length[qs]}],{AbsoluteThickness[2],Darker[Green]}],PlotLegends->SwatchLegend[Append[Map["ϱ("<>ToString[#]<>",x)"&,qs],"Least squares"]]];
Show[{grData,grQReg},ImageSize->400]
AutoCollapse[]

We can see that the new regression quantiles computed for 0.05, 0.25, 0.5, 0.75, and 0.95 have not changed significantly:

Grid[List/@Flatten[Transpose[{Chop[qrFuncs],Chop[qrFuncsAlt]}]],Alignment->Left,Dividers->{{True,True},Flatten@{Table[{True,False},{Length[qrFuncs]}],True}}]
AutoCollapse[]

ant that they are still good for separating the un-altered data:

tbl=Table[{qs[[i]],Length[Select[data,#[[2]]≥(qrFuncsAlt[[i]]/.x->#[[1]])&]]/Length[data]//N},{i,Length[qs]}];TableForm[tbl,TableHeadings->{None,{"quantile","fraction\nabove"}}]
AutoCollapse[]

Also we can see that the least squares fit of (11) has significantly changed:

fFunc

fFuncAlt

Quantile regression with B-splines

We are going to consider two implementations of Quantile Regression (QR) calculation with B-spline bases. The first implementation is based on the Linear Programming (LP) formulation of the quantile minimization problem. The second implementation is a direct translation of the non-LP minimization formulation. (Mathematica’s functions LinearProgramming and Minimize are used respectively.)

data=N[SinWithUpwardTrend[1000,12,100]];
ListPlot[data,AxesLabel->{"x","y"},PlotRange->All,ImageSize->400]

qs={0.05,0.25,0.5,0.75,0.95};

qrFuncsBSplines=QuantileRegression[data,20,qs,Method->{LinearProgramming,Method->"CLP"}];

fFunc=Fit[data,Join[{1,x},Table[Sin[x/i],{i,5,10}]],x]

grData=ListPlot[data,PlotRange->All];grQReg=Plot[Evaluate[MapThread[Tooltip[#1,#2]&,{Append[Through[qrFuncsBSplines[x]],fFunc],Append[Map["ϱ("<>ToString[#]<>",x)"&,qs],"Least squares"]}]],{x,Min[data[[All,1]]],Max[data[[All,1]]]},PlotStyle->Append[Table[{AbsoluteThickness[2],Blend[{Yellow,Darker[Red]},i/Length[qs]]},{i,Length[qs]}],{AbsoluteThickness[2],Darker[Green]}],PlotLegends->SwatchLegend[Append[Map["ϱ("<>ToString[#]<>",x)"&,qs],"Least squares"]]];
Show[{grData,grQReg},ImageSize->400]

Outlier detection

Data

See the definition in MathWorld : http://mathworld.wolfram.com/Outlier.html .

normdata=Append[RandomReal[NormalDistribution[],100],-6];
g1=Histogram[normdata,PerformanceGoal->"Speed"]

g2=ListPlot[Insert[Table[i^2+1000Random[],{i,100}],6000,30],PlotStyle->{Red,AbsolutePointSize[4]},PerformanceGoal->"Speed"]

points=Cases[g2,Point[x___]:>x,∞][[1]];
Shallow[points]

Outlier detection

Finding the outlier

qFuncs=QuantileRegression[points,2,{.95}];

ListPlot[{qFuncs[[1]]/@points[[All,1]],(qFuncs[[1]]/.x->#&)/@points[[All,1]],points}]

Estimation of conditional density distributions

Set-up

Assume we have temperature data for a given location and we want to predict today’s temperature at that location using yesterday’s temperature. More generally, the problem discussed in this section can be stated as “How to estimate the conditional density of the predicted variable given a value of the conditioning covariate?”

One way to answer this question is to provide a family of regression quantiles and from them to estimate the Cumulative Distribution Function (CDF) for a given value of the covariate. In other words, given data pairs (Subscript[{(Subscript[x, i]),(Subscript[y, i])}, i=1]^n) and a value (Subscript[x, 0]) we find (Subscript[CDF, (Subscript[x, 0])])(y). From the estimated CDF we can estimate the Probability Density Function (PDF). We can go further and use Monte Carlo type of simulations with the obtained PDF’s (this will be discussed in another post).

The experiments in this section follow the example sub-sub-section “Daily Melbourne Temperatures” of the book “Quantile regression” by Roger Koenker.

Consider the temperature time series of Atlanta, GA, USA from 2006.01.01 to 2014.01.13 :

location={"Atlanta","GA"};
(*location="Melbourne";*)
tempData=WeatherData[location, "Temperature",{{2000, 1, 1}, {2014, 10, 12}, "Day"}];
tempData//Length

jan1Pos=Map[AbsoluteTime,Cases[DateList/@tempData["TimeList"][[1]],{_,1,1,___}]];
DateListPlot[tempData,Joined->True,PlotRange->All,Axes->False,Frame->True,GridLines->Automatic,FrameLabel->{None,"temperature, °C"},PlotLabel->If[ListQ[location],Row[Riffle[location,", "]],location],AspectRatio->1/3,ImageSize->800]

We can see that this time series is heteroscedastic -- the range of temperatures is wider in the winter months.

Using the time series data let us make pairs of yesterday and today data:

tempPData=Partition[tempData["Values"],2,1];
tempPData=tempPData/.Quantity[q_,_]:>q;
tempPData=Select[tempPData,VectorQ[#,NumberQ]&];
tempPData//Dimensions

(We pair up two consecutive temperatures.)

ListPlot[tempPData,PlotRange->All,Axes->False,Frame->True,GridLines->{FindDivisions[{##},10]&,FindDivisions[{##},10]&},GridLinesStyle->Directive[GrayLevel[0.8],Dashed],AspectRatio->1,FrameLabel->{"yesterday's temperature, °C","today's temperature, °C"},PlotLabel->If[ListQ[location],Row[Riffle[location,", "]],location]]
AutoCollapse[]

Estimation of conditional density distributions

Regression quantiles

Let us calculate the regression quantiles for a number of quantiles.

The considered set of quantiles is qs:={0.02,0.05,0.01,…,0.95,0.98}:

qs=Join[{0.02},FindDivisions[{0,1},20][[2;;-2]],{0.98}]//N
qs//Length

Next we calculate the regression quantiles using 3rd order B-spline basis functions over 5 knots.

AbsoluteTiming[
qFuncs=QuantileRegression[tempPData,5,qs,Method->{LinearProgramming,Method->"CLP"}];
]

gr1=ListPlot[tempPData];
gr2=Plot[Evaluate[Through[qFuncs[x]]],{x,Min[tempPData[[All,1]]],Max[tempPData[[All,1]]]},PlotPoints->130];
Show[{gr1,gr2},PlotRange->All,Axes->False,Frame->True,GridLines->{FindDivisions[{##},10]&,FindDivisions[{##},10]&},GridLinesStyle->Directive[GrayLevel[0.8],Dashed],AspectRatio->1,FrameLabel->{"yesterday's temperature, °C","today's temperature, °C"},PlotLabel->If[ListQ[location],Row[Riffle[location,", "]],location],ImageSize->400]

Given a temperature value, say (Subscript[t, 0])=2°C, we can estimate (Subscript[CDF, (Subscript[t, 0])])(t) using the quantiles and the values at (Subscript[t, 0]) of the corresponding regression quantile functions.

t0=15;
xs=Through[qFuncs[t0]];
cdfPairs=Transpose[{xs,qs}];

We can get a first order (linear) approximation by simply connecting the consecutive points of (Subscript[{(Subscript[x, i]^(Subscript[t, 0]))(Subscript[q, i])}, i=1]^(|qs|))  as it is shown on the following plot.

ListLinePlot[cdfPairs,PlotMarkers->Automatic,Axes->False,Frame->True,PlotLabel->"Estimated CDF for "<>ToString[t0]<>"°C",ImageSize->400]

Since we want to also derive (Subscript[PDF, (Subscript[t, 0])])(t) from the estimated (Subscript[CDF, (Subscript[t, 0])])(t), make plots, and do other manipulations, it is better to make and use an interpolation object for (Subscript[CDF, (Subscript[t, 0])])(t):

qCDFInt=Interpolation[cdfPairs,InterpolationOrder->1];

Here is a plot with both (Subscript[CDF, (Subscript[t, 0])])(t) and the derived from it (Subscript[PDF, (Subscript[t, 0])])(t).

Plot[{qCDFInt[x],qCDFInt'[x]},{x,Min[xs],Max[xs]},GridLines->{Prepend[Map[{#,Dashed}&,Pick[xs,MemberQ[{0.05,0.25,0.5,0.75,0.95},#]&/@qs]],{t0,GrayLevel[0.8]}],None},Axes->False,Frame->True,PlotLegends->SwatchLegend[{"CDF","PDF"}],PlotLabel->"Estimated CDF and PDF for "<>ToString[t0]<>"°C",ImageSize->400]

On the plot the dashed vertical grid lines are for the quantiles {0.05,0.25,0.5,0.75,0.95}; the solid vertical gray grid line is for (Subscript[t, 0]).

Estimation of conditional density distributions

Bulk comparison

Let us define a functions of the CDF and PDF approximation and plots:

CDFEstimate[t0_]:=CDFEstimate[qs,qFuncs,t0];
CDFEstimate[qs_,qFuncs_,t0_]:=
Interpolation[Transpose[{Through[qFuncs[t0]],qs}],InterpolationOrder->1];

Clear[CDFPDFPlot]
CDFPDFPlot[t0_?NumberQ,qCDFInt_InterpolatingFunction,qs:{_?NumericQ..}:qs,opts:OptionsPattern[]]:=
Block[{qsGL={0.05,0.25,0.5,0.75,0.95},xsGL},
xsGL=Pick[Flatten@qCDFInt["Grid"],MemberQ[qsGL,#]&/@qs];
Plot[{qCDFInt[x],qCDFInt'[x]},{x,qCDFInt["Domain"][[1,1]],qCDFInt["Domain"][[1,2]]},GridLines->{Prepend[MapThread[{#2,{Dashed,Blend[{Blue,Green,Pink},#1]}}&,{qsGL,xsGL}],{t0,GrayLevel[0.6]}],None},PlotRange->{0,1},Axes->False,Frame->True,PlotLabel->"Estimated CDF and PDF for "<>ToString[t0]<>"°C",opts]
];

Using these definitions let us plot the estimated CDF’s and PDF’s for a collection of temperatures. Note that the quantiles are given with vertical dashed grid lines; the median is colored with green. The values of the conditioning covariate are given on the plot labels and are marked with a solid, gray vertical line.

Grid@Map[CDFPDFPlot[#,CDFEstimate[#],ImageSize->200]&,{{18,20,22,24},{26,28,30,32},{34,35,36,37}},{-1}]

Grid@Map[CDFPDFPlot[#,CDFEstimate[#],ImageSize->200]&,{{-1,0,5,7},{10,12,15,17},{20,22,25,27}},{-1}]

In the grid of CDF-and-PDF plots we can see that:
1. when the temperature is lower there is higher probability that the next day the temperature is going to be higher,
2. when yesterday’s temperature (Subscript[t, 0]) is within the range [10,25]°C the median almost coincides with (Subscript[t, 0]),
3. for medium and high temperatures we have distributions similar to the skew normal distribution.

Estimation of conditional density distributions

Let us look at the results in the previous slide from a different perspective using mosaic plots.

tvals=tempData["Values"]/.Quantity[q_,_]:>q;

Make a piecewise function using a grid of quantile values:

qvals=Quantile[tvals,Range[0,1,1/11]];
qrf=QuantileReplacementFunc[qvals];

Make a piecewise function using a uniform grid

rvals=Rescale[Range[1,12],{1,12},{Min[tvals],Max[tvals]}];
rrf=QuantileReplacementFunc[rvals];

Compare the tow labeling functions:

Grid[{{"Quantiles grid","Uniform grid"},{qrf,rrf}},Dividers->All]

Compare the mosaic plots corresponding to the two labeling functions:

Grid[{{"Quantiles grid","Uniform grid"},{MosaicPlot[Map[qrf,tempPData,{-1}],ImageSize->350],MosaicPlot[Map[rrf,tempPData,{-1}],ImageSize->350]}},Dividers->All]

Combination with sequence mining techniques

Consider the temperature of Atlanta

location={"Atlanta","GA"};
(*location="Melbourne";*)
tempData=WeatherData[location, "Temperature",{{2010, 1, 1}, {2014, 10, 12}, "Day"}];
tempData//Length

Make a data array:

tempDataArray=Transpose[{tempData["TimeList"][[1]],tempData["Values"]/.Quantity[x_,_]:>x}];

Time series plot:

grDLP=DateListPlot[tempDataArray,Joined->True,PlotRange->All,Axes->False,Frame->True,GridLines->Automatic,FrameLabel->{None,"temperature, °C"},PlotLabel->If[ListQ[location],Row[Riffle[location,", "]],location],AspectRatio->1/3,PlotRange->All]
AutoCollapse[]

Find regression quantiles:

qs=Join[{0.02},FindDivisions[{0,1},6][[2;;-2]],{0.98}]//N
qs//Length

AbsoluteTiming[
qFuncsDLP=QuantileRegression[tempDataArray,20,qs,Method->{LinearProgramming,Method->"CLP"}];
]

Plot regression quantiles and time series data:

AbsoluteTiming[grQFuncsDLP=Plot[Evaluate[Through[qFuncsDLP[x]]],{x,Min[tempDataArray[[All,1]]],Max[tempDataArray[[All,1]]]},PlotPoints->130];
]

Show[{grDLP,grQFuncsDLP}]

Combination with sequence mining techniques

Let us extract states from the fitted regression quantiles:

FindQRRange[{x_,y_}]:=
Block[{qfs,pfunc},
qfs=Through[qFuncsDLP[x]];
pfunc=QuantileReplacementFunc[qfs];
pfunc[y]
];

qstates=FindQRRange/@tempDataArray[[All]];

Here is how the state sequence looks like:

Shallow[qstates]

Here is the mosaic plot picturing the conditional probabilities of the state occurrence:

MosaicPlot[Partition[qstates,2,1]]

We can alternatively use tries:

ptrie=TrieNodeProbabilities[TrieCreate[Partition[qstates,2,1]]];

TrieForm[ptrie,AspectRatio->1/1.2]

paths=Rest/@TrieRootToLeafPaths[ptrie];
Manipulate[
TableForm[SortBy[Map[Flatten,paths],If[sortAxis≤2,{d*#[[sortAxis]],d*#[[sortAxis+2]],d*#[[sortAxis]]}]&],TableDepth->2],
{sortAxis,{1,2,3,4}},
{{d,-1,"direction"},{1->"Asc",-1->"Desc"}}]

Find fit of non-linear data

The problem

As stated in the Community site:

http://community.wolfram.com/groups/-/m/t/224105

“I am attempting to fit a function to data (pasted below) which has an s - shaped appearance. What parameter' s should I use?”

data={{0,2.52},{1.,2.83},{2.,3},{3.,3.2},{4.1,3.35},{5.,3.47},{6,3.57},{7,3.66},{8,3.76},{8.5,3.81},{9,3.85},{9.5,3.89},{10.1,3.94},{10.5,3.98},{11,4.01},{11.5,4.06},{12,4.09},{12.5,4.15},{13,4.19},{13.5,4.25},{14,4.3},{14.5,4.35},{15,4.41},{15.6,4.47},{16,4.53},{16.5,4.6},{17.,4.68},{17.5,4.77},{18,4.85},{18.5,4.96},{19,5.11},{19.55,5.34},{19.7,5.44},{19.9,5.58},{20.1,5.91},{20.3,6.27},{20.5,7.14},{20.6,7.14},{20.8,7.81},{20.9,8.32},{21,7.75},{21.2,9.07},{21.4,9.49},{21.5,9.71},{21.6,9.83},{21.8,10},{22,10.18},{22.1,10.21},{22.2,10.25},{22.3,10.27},{22.5,10.3},{22.7,10.42},{22.9,10.47},{23.1,10.52},{23.3,10.59},{23.5,10.63},{23.7,10.67},{24,10.74},{24.2,10.78},{24.4,10.8},{24.6,10.82},{24.8,10.84},{25,10.87}};

dGr=ListPlot[data,PlotRange->All,PlotStyle->Red,Frame->True]

Find fit of non-linear data

Solution

I would like to propose a different way of doing the non-linear fitting using Quantile regression with B-splines. The advantages of the approach are that it does not need guessing of the fit functions and it requires less experimentation.

After several experiments with the number of knots we can find a regression quantile using 3d order B-splines of 18 knots that approximates the data well.

qfunc=Simplify[QuantileRegression[data,18,{0.5},InterpolationOrder->2][[1]]];

qfGr=Plot[qfunc[x],{x,Min[data[[All,1]]],Max[data[[All,1]]]},PlotPoints->700];
Show[{dGr,qfGr},Frame->True]

Here is the relative error estimate:

ListPlot[Map[{#[[1]],(qfunc[#1[[1]]] - #1[[2]])/#1[[2]]}&,data],Filling->Axis,Frame->True]

Find fit of non-linear data

Possible aliasing errors

I would like to point out that if we put large number of knots we can get aliasing errors effect:

qfunc=Simplify[QuantileRegression[data,40,{0.5},InterpolationOrder->2][[1]]];

qfGr=Plot[qfunc[x],{x,Min[data[[All,1]]],Max[data[[All,1]]]},PlotPoints->700];
Show[{dGr,qfGr},Frame->True]

Find fit of non-linear data

Non-uniform grid of knots

If we select a non-uniform grid of knots according to gradients of the data we get a good approximation with one attempt.

knots=
Join[
Rescale[Range[0,1,0.2],{0,1},{Min[data[[All,1]]],Max[data[[All,1]]]3/4}],Rescale[Range[0,1,0.1],{0,1},{Max[data[[All,1]]]3/4,Max[data[[All,1]]]}]
]

qfunc=QuantileRegression[data,knots,{0.5}][[1]];

The grid lines in the plot below are made with the knots.

qfGr=Plot[qfunc[x],{x,Min[data[[All,1]]],Max[data[[All,1]]]}];
Show[{dGr,qfGr},GridLines->{knots,None},GridLinesStyle->Directive[GrayLevel[0.8],Dashed],Frame->True]

Here are the relative errors:

ListPlot[Map[{#[[1]],(qfunc[#1[[1]]] - #1[[2]])/#1[[2]]}&,data],Filling->Axis,Frame->True]

Here is the approximation function after using PiecewiseExpand and Simplify:

qfunc=Evaluate[Simplify[PiecewiseExpand[qfunc[[1]]]]]&;
qfunc

qfunc//InputForm

Future plans

Naive Bayesian classifier using Quantile regression

2D Quantile regression envelopes

More documented examples and applications

Directional quantile envelopes 2D

Problem

Given a cloud of points (2D and higher dimensions) how can we find the envelope that separates them into two sets according to a quantile parameter q?

Currently worked on

I.e. it is work in progress …

Directional quantile envelopes 2D

Let us illustrate the naive approach of finding the quantile curves in a set of uniformly angled directions (without considering Tukey depths etc.).

Two envelopes

Where to find the package and documentation

Package and documentation

[1] Anton Antonov, Quantile regression Mathematica package, source code at GitHub, https://github.com/antononcube/MathematicaForPrediction, package QuantileRegression.m, (2013).

[2] Anton Antonov, Quantile regression through linear programming, usage guide at GitHub, https://github.com/antononcube/MathematicaForPrediction, in the directory “Documentation”, (2013).

Blog posts

[3] Anton Antonov, Quantile regression through linear programming, “Mathematica for prediction algorithms” blog at WordPress.com, 12/16/2013, http://mathematicaforprediction.wordpress.com/2013/12/16/quantile-regression-through-linear-programming/

[4] Anton Antonov, Quantile regression robustness, “Mathematica for prediction algorithms” blog at WordPress.com, 12/23/2013, http://mathematicaforprediction.wordpress.com/2013/12/23/quantile-regression-robustness/

[5] Anton Antonov, Quantile regression with B-splines, “Mathematica for prediction algorithms” blog at WordPress.com, 1/1/2014,  http://mathematicaforprediction.wordpress.com/2014/01/01/quantile-regression-with-b-splines/

[6] Anton Antonov, OPL code for quantile regression with B-splines, “Mathematica for prediction algorithms” blog at WordPress.com, 1/10/2014, http://mathematicaforprediction.wordpress.com/2014/01/10/opl-code-for-quantile-regression-with-b-splines/

[7] Anton Antonov, Estimation of conditional density distributions, “Mathematica for prediction algorithms” blog at WordPress.com, 1/13/2014, http://mathematicaforprediction.wordpress.com/2014/01/13/estimation-of-conditional-density-distributions/

[8] Anton Antonov, Find Fit for Non-linear data,  “Mathematica for prediction algorithms” blog at WordPress.com, 4/19/2014, http://mathematicaforprediction.wordpress.com/2014/04/19/find-fit-for-non-linear-data/

Other

[9] Anton Antonov, Mosaic plots for data visualization, usage guide at GitHub, https://github.com/antononcube/MathematicaForPrediction, in the directory “Documentation”, (2014).

[10] Anton Antonov, Tries, usage guide at GitHub, https://github.com/antononcube/MathematicaForPrediction, in the directory “Documentation”, (2014).

[11] Anton Antonov, Mosaic plots for data visualization,  “Mathematica for prediction algorithms” blog at WordPress.com, 3/17/2014, http://mathematicaforprediction.wordpress.com/2014/03/17/mosaic-plots-for-data-visualization/

[12] Anton Antonov, Tries with frequencies for data mining,  “Mathematica for prediction algorithms” blog at WordPress.com, 12/6/2013, http://mathematicaforprediction.wordpress.com/2013/12/06/tries-with-frequencies-for-data-mining/